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We develop a theoretical framework to understand the preparation and relaxation of a metastable 
Mott insulator state within the first excited band of a ID optical lattice. The state is loaded by 
"lifting" atoms from the ground to the first excited band by means of a stimulated Raman transition. 
We determine the efi^ect of pulse duration on the accuracy of the state preparation for the case of 
a Gaussian pulse shape. Relaxation of the prepared state occurs in two major stages: double- 
occupied sites occurring due to quantum fiuctuations initially lead to interband transitions followed 
by a spreading of particles in the trap and thermalization. We find the characteristic relaxation 
times at the earliest stage and at asymptotically long times approaching equilibrium. Our theory 
is applicable to recent experiments performed with ID optical lattices [T. MuUer, S. Foiling, A. 
Widera, and I. Bloch, Phys. Rev. Lett. 99, 200405 (2007)]. 
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I. INTRODUCTION 

Cold atoms loaded into an optical lattice are a natu- 
ral system to study many-body physics far from equi- 
librium Experiments at the University of Mainz 
and NIST have succeeded in promoting a large fraction 
of bosons from the ground band to excited lattice or- 
bitals 0, Many-body phenomena within these ex- 
cited orbitals can be strikingly different from the ground 
band; theoretical proposals include the existence of su- 
persolid, algebraic bond liquid, and p-wave superfluid 
phases [1 SH, 0, i, H El, lH ■ The inherent metasta- 
bility of the excited band is the main experimental barrier 
to realizing these models in optical lattice systems. Isac- 
sson et al. pointed out that the absence of phase space 
for initial decay can lead to long excited band lifetimes 
1/r when compared to the nearest neighbor hopping time 
[5|. Miiller et al. confirmed this prediction experimen- 
tally In a related system, Spielman et al. studied 
the preparation and relaxation of bosons promoted to 
the first excited band in an array of quasi-2D lattices Q . 
Recently, Stojanovic et al. calculated the excited band 
lifetime for a superfluid of bosons in a shallow 2D double- 
well optical lattice and found that the lifetimes could be 
thousands of time longer than nearest neighbor tunneling 
times 13] . 

As a step towards a complete theory of the metasta- 
bility of the excited band, we present a detailed char- 
acterization of the preparation and relaxation of one of 
the simplest excited band states, the ID "p-orbital in- 
sulator". We define this state and describe in detail its 
preparation for the case of a Gaussian pulse in sections 
im and mil We compute the initial decay rate F using 
Fermi's Golden Rule for a wide range of well depths in 
section IIVI On long time scales our system will evolve 
to a thermal Bose-Einstein distribution; we compute its 
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parameters in section |Vl In section IVll we study the late 
stages of relaxation to this equilibrium using a Boltzmann 
equation formalism. 

Throughout this work we use the recoil energy Ej^ = 
h^/2mX^, where m and A are respectively the atomic 
mass and wavelength of the laser creating the lattice. 
Note the optical lattice separation is A/2, since the lat- 
tice potential is proportional to the intensity of the laser. 
When we make comparisons to experiment, we use the 
parameters from 2] where Ej^ = h x 3.273 kHz. 



II. QUALITATIVE CONSIDERATIONS AND 
MAIN RESULTS 

A. Adiabatic loading of the excited band 

A condensate loaded into the ground band s orbitals of 
an optical lattice at a density of one boson per site will 
form a Mott insulator for sufficiently deep wells where 
the on-site interaction is much stronger than the tun- 
neling p^ . [isj . Starting from that ground state, every 
atom can be excited to the orbital of the first excited 
band. In a typical optical lattice setup, the lattice poten- 
tial in the y— and z— directions is much stronger than 
in the x— direction; this lifts the degeneracy among the 
three p orbitals and makes it possible to freeze out exci- 
tations to the Py and pz orbitals. On time scales short 
compared to the excited band lifetime 1/F, the system is 
well described by an effective excited band Bose Hubbard 
Hamiltonian. If the nearest neighbor hopping energy ti 
is sufficiently small compared to the on-site interaction 
U, the system will have a ground state with an excitation 
gap of order U. A schematic of the p-orbital insulator is 
shown in Fig. [T] 

To be more explicit, a laser pulse causing a stimulated 
Raman transition can be used to promote bosons from 
the ground band to the p-orbital insulator [16] . The op- 
tical lattice potential is anharmonic, so the energy dif- 
ference between the ground and first excited bands eio 
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FIG. 1: (Color online.) A sketch of the p-orbital insulator at 
a density of one particle per site. For the p-orbital insulator 
the hopping energy in the excited band ti is much smaller 
than on the on-site interaction Uu . The density of virtually 
double occupied sites p « 2ff/!7ii; see Eq. (|III.2Hl . The 
darker-green orbitals indicate the dominant on-site density, 
while the lighter-green orbitals indicate virtual hops to the 
left and right. 



is larger than the energy difference between the first and 
second excited bands €21- The applied laser pulse can be 
considered adiabatic if its corresponding Rabi frequency 
is small compared to the anharmonicity eio — £21- Allow- 
ing hops of atoms between the wells of the optical lattice 
modifies the spectrum of the many-body atomic system 
by introducing the small hopping energy scale ti . 

A simple system of two atoms on two sites can eluci- 
date many aspects of that modification; we discuss this 
case in Section fill Al First, virtual delocalization of the 
atoms from their respective sites results in the transition 
frequency shifting from its nominal value by 0{t\/U). 
Second, the evolution of the system under the pulse is 
characterized by two time scales : the inverse of the on- 
site interaction 1/U and the inverse of the frequency split- 
ting U /t\. In the case of a simple two-atom system, the 
amplitude of the p-state formed by a 7r-pulse of length 
T initially grows with r reaching a local maximum at 
T '-^ jj log J- . The infidelity for a pulse with that optimal 

pulse length is of order jjt log . Further increase of 
the pulse length first increases the infidelity because of 
the optical nutation between the two-atomic p-state and 
a state in which only one of the atoms is in that state. In 
the case of two atoms the full fidelity is reached asymp- 
totically for TT-pulse lengths exceeding significantly the 
t\/U scale. 

In a lattice of many atoms, the p-state acquires a finite 
lifetime 1/F. For a typical experimental setup, T ^ J7, 
but F ~ tl/U. The latter relation makes it practically 
hard to get arbitrarily close to the ground-state of the 
p-state Mott insulator. We argue that the smallest infi- 
delity one may reach before the relaxation of the p-state 

sets in is of the order of ^ log^ ^ (see Section UlI Bp . 



B. Stages of Relaxation 

The p-orbital insulator ground state has a fraction of 
double-occupied sites of 0{ti/U)'^. The on-site interac- 
tion causes pairs of particles to make transitions out of 
the p-band, releasing large amounts (^ eio) of energy. 
Metastability is guaranteed if however the fraction of 
these doubly occupied sites is sufficiently small. In sec- 
tion llVi we present a detailed calculation for the initial 
decay of the p orbital insulator. 

We can estimate the initial decay rate F by applying 
the Fermi Golden Rule to transitions transferring two 
atoms from the p-band (n = 1) into states in band n = 2 
and in band n = 0. The rate F is proportional to the 
product of the probability of double occupancy of a state, 

(ti/U)'^, the square of the inter-band transition matrix 
element ~ C/^, and the combined density of final states. 
The latter one is ^ l/t2 and is essentially controlled by 
the wider band n — 2. The leading order result for the 
initial decay rate is 

T^tyt2. (Ill) 

Note in Eq. pi.ip that the decay rate F does not depend 
on the interaction strength to leading order. For V > 
26i?fl, the bandwidth of the second excited band becomes 
too narrow for this two-body decay process to conserve 
energy; our analysis predicts a jump discontinuity in the 
decay rate at the threshold value of well depth. A more 
complicated three-body process is the dominant initial 
decay channel for well depths beyond this threshold. 

Assuming that the atoms are thermally isolated in the 
trap, the very large energy deposited in the system by the 
initial 7r-pulse eventually is transformed into the thermal 
energy of the atoms in the trap (see section|V]for details). 
At equilibrium, the effective temperature T is on the scale 
of the highly energetic eio. Because of the high thermal 
energy per atom, the equilibrium density of the system 
is considerably lower than the initial density due to the 
spreading of atoms within the magnetic trap. We treat 
the asymptotic approach to this final equilibrium using 
a Boltzmann equation formalism in section IVll 



III. ADIABATIC LOADING OF THE EXCITED 
MANY-BODY STATE 

We want the 7r-pulse to excite N bosons from the true 
ground state to the ground state within the Px band with 
high fidelity. In order to elucidate the effect of the length 
of the pulse we first investigate a model of two sites with 
two energy levels mixed by tunneling. We then comment 
on the full iV-site problem, and show that the condition 
for adiabaticity remains Ut ^ H, in the sense that the 
density of excitations out of the px band ground state is 
small in this limit. 
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A. Two-site model 

Consider a model of two bosons in a well of two sites. 
Each site has an s level and a p level, but we only in- 
clude tunneling in the p-band. We denote the band 
splitting by eio, the excited band hopping energy by 
ti, and the on-site interaction between bosons in bands 
ni,n2 as U„-^n2- We will describe our model in the sec- 
ond quantized formalism, using the bosonic operators 
bw, bLi, fcflo, bj^j where {L, R} refers to the left and right 
positions and {0, 1} refers to the band index. Taking the 
s orbital energy to be zero, the Hamiltonian H is given 
by the expression: 



+ {L^R). 



(III.l) 



As our initial state was symmetric with respect to left- 
right exchange, the relevant Hilbert space consists only of 
parity-symmetric states. We have sketched these states 
in Fig. [21 Neglecting terms of 0{t\/U^), the six parity- 
symmetric eigenstates of H are 
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The corresponding eigenvalues for these states to the 
same order in ti/U are 
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TABLE I: A conversion table between the dimensionless pa- 
rameters U = UuT and k defined in Eq. (|III.6|l and the ex- 
perimentally relevant values of well depth V and pulse time t. 
For fixed pulse time, ft is exponentially sensitive to changes 
in well depth whereas U is nearly constant over a wide range 
of well depth energies; this difference in sensitivity allows an 
experimentalist to tune the two parameters nearly indepen- 
dently. 



with the Gaussian envelope A(i) = ^ exp(— i^/r^). We 
will choose the dimensionless pulse strength a to optimize 
a single 7r-pulse . As the excitation can be considered a 
standard two-photon process, the carrier frequency of the 
pulse Lo is chosen to be half the energy difference of states 
|1) and 1 5): 



^ — 2 (^55 — Hix 
2tj 



Uu 



(in.5a) 
(IIF5b) 



To simplify notation in the following discussion, we define 
the dimensionless time u — t/r, the dimensionless on-site 
interaction in the p-band U = Uht and the dimensionless 
detuning k between the energy of intermediate state |3) 
and the drive frequency: 
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Within our model, we want to prepare the system in 
state |1) and, by applying a 7r-pulse, adiabatically evolve 
it into state |5). This process is the two-site equivalent 
of promoting a ground band Mott Insulator into the p- 
orbital insulator. We consider a general pulse term of the 
form: 



pulse 



Ait) e 



bl^bLo + b'ji.bHo] +h-c., (in.4) 



Experimentally, k and U can be controlled nearly inde- 
pendently since k depends much more sensitively on well 
depth than U. In our calculation, we will treat these two 
parameters as independent. A conversion table between 
{k, U} and the more experimentally relevant variables of 
well depth V and pulse time r are given in Table HI 

To discuss the state's time evolution, we will write a 
Schrodinger equation for the amplitudes of each of the 
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FIG. 2: (Color online.) The six parity-symmetric basis states 
of the two-site model in section UlIAl When the excited band 
hopping energy ti = 0, these states are true eigenstates. For 
ti 7^ 0, states |3) and |4) are mixed, as are states |5) and |6). 
The eigenstates are given by Eq. (|III.2p to leading order in 
ti/U. 



six parity-symmetric eigenvectors: 

\m)=ai\l) 



as |3) + e'"'a4 |4) 



02 |2> + 

+e'^'^'a5\5) + e'^'^''aG\6) . (III.7) 



The Hilbert space can be considered as two independent 
subspaces {|1) , |3) , |5)} and {|2) , |4) , |6)} weakly cou- 
pled by the pulse. The amplitudes of the second sub- 
space are much smaller because their corresponding en- 
ergies are shifted by order the Mott gap compared to the 
energies of the first subspace. The strategy for solution 
in the k <C 1 limit consists of restricting the evolution to 
the {|1) , |3) , |5)} subspace, and then using this solution 
as a drive term to find the evolution of the much smaller 
amplitudes in the {|2) , |4) , |6)} subspace. The details 
are located in Appendix |^ 



The important parameters to determine from our 
analysis are the dimensionless strength of the optimal 
pulse a, the final amplitude of the p-orbital insulator 
|a5(cxD)|2 subject to this optimal pulse, and the ampli- 
tudes {02, 04, ae}. In the k <C 1 and U ^ 1 limit, we can 
expand these parameters in the small parameters k and 
l/il: 
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with r — Uii/Uio and the dimensionless interaction en- 
ergies {1/2,114, Ue} defined by 



U2 



■U, 



2-r 



(Ill.lOa) 
(Ill.lOb) 
(III.lOc) 



Among the exponentially suppressed amplitudes 
{a2,a4,ag}, by far the largest is ag. As shown in 
Appendix lA 41 the dominance of ag is a consequence of 
Uii being smaller than either Uio or Uqq. 

An important observation about Eq. (jIII.8P is that per- 
fect fidelity is impossible for a finite n even with f/ — > 00. 
This occurs because the spectral width of the pulse is suf- 
ficient to populate the intermediate state. The Hamilto- 
nian H + Hp^isc in the limit ?7 — )■ cxd is similar to that of a 
spin-1 system in a magnetic field, where a small quadratic 
Zeeman shift leads to nonequidistant energy levels (see 
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eas lIII.ll and IIII.4|) . In principle, if one would increase 
the pulse length and went to sufficiently large values of 
K, one could again approach high fidelities with a single 
frequency. This limit is infeasible, however, due to the 
finite lifetime of the excited band. 

As shown in Fig. [H for a given well depth there is 
an optimal pulse time t which balances the decreasing 
infidelity at larger U against the increasing infidelity at 
larger k. We can approximate the optimal pulse by using 
the following rough form for 1 — | as | : 

l-\a,\' ^ae-^^- + b^T', (III.ll) 

with a,b ^ 0{1). The first term comes from the suppres- 
sion of |a6p, and the second from the growth of |aip and 
ag. The asymptotic expansion of the optimal pulse time 
Topt in Eq. pil.lip for small ti/U is 

ropt = |log^ + o(iloglog|). (III.12) 

The expression for the resulting infidelity is 

|a5(ropt)|'«|Jlog^|. (III.13) 

A two-tone pulse could make 1 05(00)!^ arbitrarily close 
to unity through the application of two 7r-pulses, but for 
small K the fidelity from a single-tone pulse is already 
better than the 80% transfer efficiency observed in ex- 
periment . The observed infidelities are more than two 
orders of magnitude larger than those calculated here. 
This striking difference could potentially come from the 
experimental pulse shape not being close to Gaussian or 
the internal structure of the Rb atom interfering with 
pulse fidelity. 

B. Full N-site Adiabatics 

We now consider the problem of loading the excited 
band in an system of N sites, modeling each site as a 
two-level system with energy splitting eio and with tun- 
neling of energy ti permitted between the excited levels 
of neighboring sites. Additionally, we make the impor- 
tant simplification of assuming all the Hubbard U pa- 
rameters are equal. Initially, the system is prepared in 
its true ground state: one particle per site, each site in 
the lower energy level. An excitation pulse A(i)e*'^* pro- 
motes bosons from the ground band to the first excited 
band (see Fig. [4]). 

We are interested in pulses sufficiently long that U ^ 
1, so that the excitation density above the ground state 
is small. We refer to such pulses as "adiabatic" , although 
this is not the conventional meaning of a slow exter- 
nal process which transfers no heat to the system fl^ . 
The three important pulse parameters are its carrier fre- 
quency Lu, its dimensionless strength a, and its duration 




FIG. 3: Probability of not finishing in the two-atom p-state 
after application of a vr-pulse (1 — [asp) plotted versus Ur at 
fixed ti/U; see Eq. (|IlL8)) . The curve minimum (correspond- 
ing to maximum fidelity) is determined by two competing 
effects. Increasing Ut exponentially suppresses |a2|^, |a4|^ 
and laej^, but also increases \aif and [asl^. The minimum 

occurs at Ur ~ log{U/ti) and has value O ^^log^ see 

Eqs. (IIII.12|I and (IIII.13|l . 




FIG. 4: (Color online.) We generalize to the case of a vr-pulse 
on A*' sites by first solving for the evolution of the system 
in the absence of tunneling, and then treat tunneling as a 
perturbation in the interaction picture. 



T. We choose the carrier frequency to maximize the oc- 
cupation of the p-orbital insulator and adjust the pulse 
strength to produce a tt— pulse, leaving the pulse dura- 
tion as the most important adjustable parameter. 

We require C/n/F ^ 1 so that the pulse time r is both 
much longer than 1/U from adiabatics considerations and 
also shorter than l/F from finite lifetime considerations. 
This condition is satisfied for a wide range of well depths 
(see Fig. [S]). 

The Hamiltonian H for this system in the frame rotat- 
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FIG. 5: Ratio of on-site interaction energy to decay rate 
versus well depth. This ratio must be large to justify a 
key assumption of this work: pulse times can both be suf- 
ficiently long to treat the system in the the adiabatic regime 
(t ^ and also respect the finite lifetime of the excited 

band (r < l/F). 



ing at the carrier frequency of the pulse oj is 



H=Ho{t) + Hi, 
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H, = (eio - c.) ^ b^h, + ^ bib,, (111.14c) 
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The initial state |'0o) is given in second-quantized nota- 
tion in the soluble limit ti = by: 
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To find the evolution of the system as a perturbation 
expansion in the small parameter ti/U, we transform to 
the interaction picture: 



z- \^i) =Hj , 
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(111.16b) 



Since we are treating all the Hubbard U terms as equal, 
the on-site interaction commutes with Hq. In the trans- 



formed basis, Hi becomes 
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B]{t) ^ cos 9b\^ +i sin ebl^, 



with the parameter given by 
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By changing to the interaction picture, we have trans- 
formed the time dependence of the pulse into time depen- 
dence of the tunnehng and detuning eio — uj. We choose 
the duration of the pulse t to correspond to a 7r-pulse 
when the frequency is resonant for a single site lu = eio- 
There are two parameters we would like to determine: 
the probability Pg that the final state after application 
of the pulse is the p-orbital insulator and the "density of 
excitations" pcxc of the final state. 

We assume for simplicity that every particle is pro- 
moted to the p-band. Since an excitation above the Ri- 
band ground state costs approximately U in energy for 
ti/U <C 1, we can define pcxc as [18] 
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(III.19) 

We note that Eq. piI.lQp does not simply count the den- 
sity of doubly excited states, since the p-orbital insulator 
has 0{tl/U'^) doubly excited states. Rather Eq. (|III.19p 
measures the relative potential energy of the final state 
in units of the on-site interaction U . 

We cannot simply apply the perturbation methods of 
the previous section to the A''-site model because the ex- 
pansion parameter changes from ti/U to ti^/N /U, which 
leads to a series that is not obviously convergent in the 
thermodynamic limit. Instead, we approach the problem 
using only intensive quantities at every stage. Consider 
turning on the tunneling between nearest neighbor sites 
one at a time. The Hamiltonian Hn with hopping allowed 
only between the first n nearest neighbor pairs is 

n 

H.n{t) E {b]+,B, + h.c.) + (eio - ^) E Bj^i 



U 



^Y.Y.^njbljbn'jbnj. 
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To compute the probability Pg that the system ends in 
^ ■ the ground state and the excitation density poxc, we de- 



rive recursion relations describing how these quantities 
change when tunneling is activated between one addi- 
tional nearest neighbor pair. Solving these recursion re- 
lations, we take n — > iV to derive results that remain 
valid in the thermodynamic limit where ti^/N/U is a 
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large parameter. This removes many of the difficulties 
present in Rayleigh-Schrodinger perturbation theory. 

We begin our discussion with |^'g), the ground state 
wave function of Hn{oo): 
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We see from Eq. (|III.21I) that the virtually double oc- 
cupied fraction is 0{t\/U'^) as illustrated in Fig. [TJ As 
discussed in Appendix [C] we can determine the ground 
state wave function by enforcing the adiabatic pulse con- 
dition: 
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where l^*") is the final many-body wave function in the 
interaction picture: 



1*7) = exp [-iT J dt'Hr^it')^ , (HI. 
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and T is the Dyson time-ordering operator. Expanding 
Eq. pil.23[) to leading order in gives 

|vi/?)^A„(oo)|vi/;') 
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As shown in Appendix [Cl the frequency u) which max- 
imizes the occupation of the p-orbital insulator is given 
by: 
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[2 cos^ d(u) cos W(y!) -f sin 2B(u) sin 2B(v!)\ . 

(IIL25c) 

We have plotted w — eio in units of t\lU in Fig. [6l In 
the limit of large U , lo eio — 2t\/U , just as in the two- 
site case considered in section UlI Al For this particular 
frequency, 




FIG. 6: The optimal pulse carrier frequency a; — eio in units of 
ti/U for the specific case of a Gaussian pulse; see Eq. (|III.25|I . 
The optimal frequency maximizes the occupation of the p- 
orbital insulator after completion of the 7r-pulse. As expected, 

for (7r 3> 1 we recover the two-site result that lj ~ eio rr- 
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(in.26a) 

2 



du I dw'e-*^"^ sin 26*6^'^"' 



(in.26b) 

U{u'-u) e{u')x 



oo pU 

du / du' sin 

-oo J — oo 



cos26'(u') -t- cos [29{u) - 29{u')] 



(III.26c) 



As expected, |5B(l)(oo)|^ \6C^^Hoo)\^ and \5D^'^^ {oo)\'^ 

are all exponentially suppressed for large U. From Ap- 
pendix [cl the recursion relations for Pg and pcxc are 



pn+l 

log - ^ 



pn 
S 



'■1 

2t|T2 

~c7^ 



<5B(i)(oo) + (5C(i)(oo) 



(in.27a) 



^2 2 
Poxc(?l + 1) =/Ocxc(n) + (5^(1) (oo) 



+2 2 

1 ,5^(1) (oo) 



(III.27b) 



Solving the recursion relations Eq. pil.27p for n ^ N 
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gives 

Pexc(^) 



Pg{N) =exp 



(55(1) (oo) ^ JC^i^oo) 

(III.28a 



+ 0' 1 1 



[/2 ' 

(III.28b) 



The exponential form of the ground state probabil- 
ity clearly is consistent with the linked cluster pic- 
ture that an extensive perturbation theory exponenti- 
ates. For the specific case of a Gaussian envelope, 
A{t) = ^expi-t/rf, the formulas for \SB'-^\oo)\^ , 
|5C(i)(oo)|^ and |i:'(i)(c5o)|^ are : 



D(i)(oo) 



due 



iUu—u 



sm 



Trerf ■ 



(IIL29a) 





/•oo 


— TT 


y ' 

J —00 




POO 


^ 2 


J —OQ 


X i 2 cos^ 




Trerf 


— cos — — 



due 



du 



iUu—u 



( Trerf 



V 2 



- (1 -I- erf u) 

4 



■ cos ■ 



(III.29b) 

C7(u' - u) 
Trerf ul 



(III.29c) 



We have plotted the sum 

Fig. [T] As shown in Appendix [Bl the excitation density 
Pcxc can be estimated for J7 3> 1 using the stationary 
phase approximation: 



4 



^ ^ 2 



with Fstat(fc) given by: 

^tat(t^) = T" 



>;tat(f/) + i;tat(-«7)* 



(III.30) 
(III.31) 



exp 



-|C7|Jl0g-^X 



1 



21og^ 4V^ 



(III.32) 



A comparison between the stationary phase estimate and 
the actual numerical integral is graphed in Fig. [71 
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FIG. 7: (Color online.) Leading order coefficient 
|5_B('^^(oo)P -|- |(5C'''^''(oo)|^ of the excitation density pcxc in 
powers of plotted versus Ut for the specific case of a 

Gaussian pulse; see Eqs. HIII.28P and pil.29|l . The stationary 
phase estimate (dashed red) of the coefficient is a very good 
approximation to the full numerical calculation (solid blue) 
for ?7r S> 1. Both curves are formally accurate only for large 
Vt. 



To conclude this section, we note that there are many 
similarities between the iV-site problem and the two-site 
problem. This is natural as we have expanded the Hamil- 
tonian to only leading order in ti/U ^ essentially linking a 
given site only with its nearest neighbors. Further neigh- 
bor terms in the Hamiltonian only appear at higher or- 
ders in ti/U, and are correspondingly suppressed. We 
similarly expect that the limitations on the pulse fidelity 
coming from a competition between U and k becomes 
important at 0{t\/U'^). 



IV. INITIAL STAGE OF DECAY 

The p-orbital insulator will decay since it is coupled 
to a continuum; nevertheless, its metastability is guaran- 
teed by the anharmonicity of the optical lattice potential 
as shown in Fig. [51 For typical experimental parame- 
ters, the dominant decay mechanism roughly consists of 
an Auger process where two bosons in the n = \ band 
decay into a nearly free particle in the n — 2 band and a 
localized n = particle at the site of interaction. We cal- 
culate this initial decay rate using Fermi's Golden Rule. 
Throughout this section, we are interested only in suffi- 
ciently deep well depths that ti/U ^ 1 and the width 
of the lowest band can be completely neglected. See Fig. 
[HI for a plot of relevant ratios of the hopping energy in 
the ground and first excited bands io, ti to the on-site in- 
teraction energies Uqq,Uiq and Un. For lattices deeper 
than 26-Efl, the width of the n — 2 band is too narrow to 
compensate for the anharmonicity of the potential and 



9 




FIG. 8: (Color online.) The principal on-site decay process in 
relatively shallow lattices {V < 26iJij) consists of two bosons 
in the n = 1 band decaying into a localized n — particle and 
an a delocalized n = 2 particle, leaving a localized hole in the 
n = 1 band. For deep lattices, the decay process is prevented 
by the anharmonicity of the energy levels. 




10 20 30 40 
Well Depth (units of E^) 



60 



FIG. 9: (Color online.) Ratios of the nearest neighbor hop- 
ping energy in bands n = 0, 1 to relevant on-site interaction 
energies plotted versus well depth. 



this simple Auger process is energetically forbidden. In- 
stead, a more complicated three-body process mediated 
by a virtual intermediate state becomes important. 



drive transitions between bands: 

^^intra = X! X! ^^^"i + X! X! ^lij^nj-l + h.C. 
n j n j 

+ I E E bUbl.b^A^,, (IV.la) 



E E E ^'^mjbnj-l + h.C. 



+ ^E E f^™«E^l^l^™^-^"^' (iv.ib) 

with tmn being the hopping matrix element from band 
m to band n and U^^^ being the interaction matrix el- 
ement which drives pairs of bosons from bands k,l to 
bands m, n. We have neglected further neighbor hop- 
pings, which are small for the bands relevant to our cal- 
culation. To simplify notation, we use shorthand forms 
for the diagonal terms i„ = t„„ and Umn = C^m"- The 
off-diagonal tunneling terms tmn only appear in pertur- 
bation theory with energy denominators 0{€m — for 
the lattice parameters we are considering these are small 
and can be neglected. The system is in a sufficiently 
dilute limit that we can restrict our analysis to s-wave 
scattering a^/A <C 1 [19]. The matrix elements J7^'„ are 
given by the expression: 



Ul,=9l Wt{z)Wt{z)Wra{z)Wn{z)dz, (IV.2a) 



^Vtrans^^>)4^^ , (IV.2b) 



where A is the wavelength of the optical lattice lasers, Wm 
are the Wannier functions for the lattice potential, and 
^vtrans jg grouud State Wannier function associated 
with the tightly confined transverse directions. We will 
work using the scaled position variables z = and 
scaled momentum variables k = in these variables, 
the first Brillouin zone covers the range — 1 < A: < 1. 

The nearest-neighbor hopping energy tm and interac- 
tion matrix element ?7jJ" obey the following scaling re- 
lations: 



J jmn 

Er 



=fn 



V 

V_ Os 
Er' a 



(IV.3a) 
(IV.3b) 



Multi-band Bose-Hubbard Hamiltonian 



Bosons loaded into an optical lattice are approxi- 
mately described by the multi-band Bose Hubbard model 
(MBH) 0]. As we want to simulate the dynamics of 
bosons promoted to the p orbital, it is useful to divide the 
MBH Hamiltonian into intraband terms TJjntra that pre- 
serve band occupation and interband terms i/intcr which 



where /„„ and cjmn are system-independent functions. 
For the experimental parameters in [2], a^/A = 0.069 . . .; 
we use this value for all subsequent calculations. 

For the following discussion, we introduce second quan- 
tized operators 6]^^, which create a localized particle in 

band n at site j, and fej^^,, which create a Bloch wave state 
in band n with quasimomentum k. We exclusively use 
i, j' for site indices, and fc, q, k' , q' for momentum indices. 
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B. Energetics 

To compute the short time decay rate of the p-orbital 
insulator, we calculate interband transitions between 
eigenstates of iyintra using Fermi's Golden Rule; see 
Eq. pV.ip . To leading order in ti/U, the p-orbital j'f/'o) 
insulator state is given in second quantized notation by: 



1 - 



-E {blb^,,+bl,bu') 



(IV.4) 

where {i'j') indicates a sum over nearest neighbors. For 
experimentally relevant parameters, the p-orbital insu- 
lator can directly decay into only two energetically per- 
missible final states, which we denote {ipAf) ^^nd \tpBf)- 
In terms of the second quantized operators bl^j and 6^^^,, 
these states are approximately given by: 

k, q)) « bl^bl^bikbij |?/'o) , (IV.5a) 

l^BfU, fc, k', q)) « bl^blb^kbik' Ho) . (IV.5b) 

The difference between {ipAf) and IV'b/) is the nature 
of the holes in the n — I band: I'lpAf) has an n = 1 
hole localized around the site of the interband transi- 
tion, whereas IV's/) has two delocalized holes. A de- 
tailed calculation gives that Decay Channel 1, where our 
initial state decays into li'Af), is allowed energetically for 
< V < 25.97 and Decay Channel 2, where our initial 
state decays into is allowed for < F < 30.93. 

The rates vanish above these threshold values of well 
depth. 

In Appendix [Dl we demonstrate how (?As/|i?intor|'*/'o) 
is suppressed by OitijU) compared to ('0yi/|^^intor|V'o) 
by solving a representative model. This result is a conse- 
quence of the n = 1 hole in IV'A/) being localized around 
the decay site j = 0. The state has a small 0(ti/C/) 

overlap with the decay site since it is orthogonal to IV'a/)- 
This difference in the overlap between final states and the 
decay site leads to the difference in the size of the matrix 
elements. We neglect i?2, the decay rate corresponding 
to Decay Channel 2, since it is Oiti/UY smaller than 
the decay rate corresponding to Decay Channel 1. 
To simplify notation, we simply refer to i?i as the initial 
decay rate F. 

C. Fermi's Golden Rule Calculation 

We now compute the principal decay rate Ri using 
the full Hamiltonian in Eq. pV.ip . The final state is 
\^Af{j, k, q)) defined in Eq. (|IV.5[) . The transition matrix 
element Ufi depends on both the quasimomentum of the 
n = 1 band hole k and the quasimomentum of the n = 2 
band particle q: 

/oo 
dz e'(''-'i>U2q{z)u\k{z)W^{z)W^{z), 
-OC 

(IV.6) 



where Unk are the Bloch functions and Wn{x) are the 
associated Wannier functions The matrix element 

between the initial and final state is 



{'^f{3,k,q) |iJi„ter| V'o) 



-AtlUflik, g)e*'^J('=-«) cos (fcTr) 



(IV.7) 



We treat the interaction between the delocalized n = 
2 boson and the p-orbital insulator using a mean field 
approximation. This gives an additional mean field shift 
to the dispersion in the n — 2 band: 

£2, -e2q + C/i2(g), (IV.8a) 

/OO 
dz\u2k{z)\' Wl{z-jnf . 

-°° j=-oo 

(IV.8b) 

The dependence of Ui2{k) on k is weak and can be ap- 
proximated by Ul2- Substituting these terms into the 
Fermi Golden Rule formula for the initial decay rate gives 



T = —tl / dkcoa^ink) dq '^^ 
Jo Jo 



eife). (IV.9) 



To get insight into Eq. (|IV.9|) . we restrict band hop- 
pings only to nearest neighbors, a valid limit for deep lat- 
tices. This approximation makes it possible to evaluate 
the decay rate F analytically. As we show here, the rate 
predicts a jump discontinuity to zero above the threshold 
well depth. Keeping only nearest-neighbor hoppings, the 
dispersions in bands n = 1, 2 are given by: 



~ h- 2i2 cos(7rfc), 
eifc ~ ei + 2ti cos(7rfc). 



(IV.lOa) 
(IV. 10b) 



Note e2k includes the mean field shift. We treat the n = 1 
band hopping energy ii as much smaller than the effec- 
tive n = 2 band hopping energy t2; this is consistent 
with experimental values. Using the simplified band dis- 
persions from Eq. (jIV.10[) . the expression for the short 
time decay rate F is 



F = 



'W7 



dk- 



U^f{k,q)\ cos^ (irk) 



" (Uiif Jl- {A- B COS nkf 



6(1- \A- Bcosnk]) , (IV.ll) 

with the parameters {A, B} only determined by well 
depth: 



A{V) = 



£2 



eo - 2ei 



2^2 



B{V) ^h/t2 



(IV.12a) 
(IV. 12b) 
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Well Depth (units of E^) Well Depth (units of E„) 



FIG. 10: (Color online.) Comparison of the model short time 
rate with a full immerical calculation of the rate for the p- 
orbital insulator to decay to \ipAf)', see Eqs. (|IV.5|l . (|IV.9|) and 
PV.III) . The model captures the presence of a logarithmic 
divergence in the rate due to a Van Hove singularity, as well 
as a jump discontinuity at threshold. 



The rate T vanishes above a threshold value 14, deter- 
mined by the condition |A(Vt) — B{Vt)\ = 1. Below 
threshold, the rate can be approximated by Eq. (jll.l[) 
except near the well depth T^ing that satisfies |A(V^ing) + 
B{Vsing)\ = 1- Here the rate diverges logarithmically due 
to the presence of a van Hove singularity. At threshold, 
the rate has a jump discontinuity to zero: 



T{V+) =0. 



(IV.lSa) 
(IV. 13b) 



The divergence in the rate and the jump discontinuity at 
threshold are also evident in the numerical evaluation of 
Eq. (jIV.9P : see Figs. [TUl and [TT] for a comparison with 
the simplified model. A full discussion of the two-body 
decay rate near threshold is given in Appendix [E] 



D. Decay of p-orbital Insulator above threshold 

Above threshold, the decay process to l^/^yi/) becomes 
energetically forbidden and one must go to higher order 
processes to determine the initial decay rate. Neglecting 
the effects of the magnetic trap, the dominant short time 
decay process is an effective three-body decay rate: 



<5 (e4fc, + 2eo - 2ei - euO , (IV. 14) 



where e^ki includes the mean field shift. This rate is en- 
ergetically allowed for V ^ iQE^. Although the second- 
order perturbation theory expression for the three-body 



FIG. 11: (Color online.) The short time lifetime in units 
of the hopping time thop for both the full and model systems; 
see Eqs. (|rV9)l and PV.11|I and Fig. [lOl 



interaction matrix element includes a sum on all inter- 
mediate bands, the contribution from the n = 2 band is 
most important: 



, -N — 8ti cos(7rfci) 



_i 2 ei -f eifei — e2fe — Cq 



(IV. 15) 



A recent paper addressed an effective three-body decay 
mechanism utilizing an intermediate virtual state in a 
related system [2^ . 



E. Combining two- and three-body decay rates 

We expand the behavior of both the two- and three- 
body short time decay rates near the threshold value Vt in 
Appendix IeI For V <Vt, we can approximately describe 
the early decay dynamics by just the two-body decay 
rate; see Eq. (IIV.9p . As y — > Vf, we find that the two- 
body rate rises linearly to a constant before jumping to 
zero: 



7 dg{V) 



16 dV 



Vt 



9(V) =^ (e~2 



2t-2 



60 - 2ei 



2U 



, (IV. 16a) 
, (IV. 16b) 



where V is the well depth in units of recoil energy. For 
well depths above the two-body threshold {V > Vt) we 
can describe the dynamics using the effective three-body 
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FIG. 12: (Color online.) Comparison of the short time three- 
body decay rate 1 + 1 + 1 + + 4 and the short time 
two-body decay rate 1 + 1 -> 0+2; see Eqs. (|IV.9|) and pV.14|l . 
Note the two-body rate exhibits a discontinuous jump to zero 
at threshold, whereas the three-body rate diverges. 

decay rate; see Eq. (jIV.14[) . As y Vf^ , the three-body 
rate becomes singular as one approaches threshold: 



r « 



8(ti 


yi/2 




2 






dk 




d 

dV 


g{v) 


(I 



This singularity can be understood in the following way: 
the virtual intermediate state is progressively longer lived 
for well depths closer and closer to threshold. At thresh- 
old the state becomes truly long-lived and there is a di- 
vergence in the decay rate. Comparisons of the short 
time two-body decay rate 1 + 1^0 + 2 and the three- 
body decay rate 1 + 1 + 1^0 + + 4 are plotted in Figs. 
[Ulandlini 

This approach only fails for a small region near thresh- 
old, where the analysis is complicated by the presence 
of an allowed real intermediate state for V < Vt- To 
benchmark the width of this inapplicable region, we cal- 
culated the well depth V;3bod where the three-body decay 
rate equals the two-body rate just below threshold (see 
Fig. [111). We found Fabod - Vt = Q.HEr, a small value 
compared to the experimentally accessible range of well 
depths. 

Linking the short time rates to the measurements made 
by Miiller et al. is made difficult because our approach 
neglects the spreading of atoms in the trap and the effect 
of an imperfect 7r-pulse fl6|. We can, however, define a 
rough time scale for the experiments where the optical 
lattice system crossed over from early stage dynamics to 
late stage dynamics. Fitting the experimental data with 
a double exponential, we can define the inverse of the 
faster decay rate to be the "short time decay lifetime". 
For the well depths considered, this lifetime was typically 



FIG. 13: (Color online.) Comparison of the short time three- 
body and two-body decay lifetimes in units of the hopping 
time thop = 2ti/h; see Eqs. (irv\9|l and (|IV.14|I and Fig. [H 
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FIG. 14: (Color online.) Comparison of the short time three- 
body and two-body decay rates near the threshold value Vt; 
see Eqs. (IIV.16al ) and (IIV.17[) . Just below threshold, the two- 
body decay rate reaches a value of about Q.QlQER/h, after 
which it jumps to zero. The three-body rate reaches this 
same value at a well depth only 0.331?^ above threshold. We 
expect that the Golden Rule approximation breaks down only 
in a small region around threshold, with a width of about Er 
on either side. 



tens of milliseconds. In Fig. [TSlwe compare this observed 
lifetime with the Fermi Golden Rule rate for the p-orbital 
insulator. Our predicted short time rates are of the same 
order as those observed in experiment for well depths 

V — 10i?fl and 20i?fl, but they differ substantially for 

V = ZQEr. 

One possible explanation for this difference with ex- 
perimental data is that the system has a small fraction 
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Well Depth (units of E^) Well Depth (units of E^) 



FIG. 15: (Color online.) Comparison of our theoretical short 
times rate with lifetimes extracted from experimental mea- 
surements published by Miiller et al. [l^. We have plotted 
the lifetime for the p-orbital insulator, as well as "imperfect" 
excited states with an additional 0.5%, 1.0% and 2.0% fraction 
of double occupied sites (DOS). The measured values appear 
consistent with roughly 1% DOS. Sources of these excitation 
infidelities could come from imperfect 7r-pulses. 



of double occupied sites due to an imperfect preparation. 
Neglecting interference terms, the Golden Rule rate is 
roughly the fraction of double occupied sites multiplied 
by the on-site decay rate of a single site. As shown 
in Fig. [lH the observed data seems to be consistent 
with about 1% of sites doubly occupied due to prepa- 
ration imperfections. Since the p-orbital insulator alone 
has roughly 2t\/U'li double occupied sites due to virtual 
hops, the short time decay rate of the experimental sys- 
tem will substantially differ from theory only when the 
fraction of "accidentally" double-occupied sites becomes 
comparable to 2tf/Ufi. As shown in Fig. [161 the ratio 
2i?/C/ii > 0.01 for V < 22Er, but rapidly falls below 
0.01 for deeper well depths. Reducing the number of 
double occupied sites through improved 7r-pulses could 
substantially improve excited band lifetimes, particularly 
for well depths V > 20Er. 



V. THERMAL EQUILIBRIUM 

After the initial stage of decay, the p-orbital insula- 
tor thermalizes and spreads in the trap. This stage of 
decay is difficult to model. At very long times the distri- 
bution is described by Bose-Einstein statistics, with the 
chemical potential fi and the temperature T'^ determined 
by total energy and particle number. The interactions 
among the particles (of order U) are much smaller than 
the temperature (of order the band splitting eio); see Fig. 

The magnetic trap reduces the equilibrium tempera- 



FIG. 16: (Color online. )The p-orbital insulator has roughly 

a fraction of double occupied sites. For V < 20Er, this 

^11 

fraction is larger than the estimated fraction of sites doubly 
occupied due to pulse infidelities in experiment 0], and we 
expect that the short time decay rate of the experimental 
system will be similar to that of the p-orbital insulator. For 
deeper well depths V > 20Eii, we expect that the short time 
decay rate is dominated by decays of sites doubly occupied 
due to pulse infidelities. 
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FIG. 17: (Color online.) Plot of on-site interaction energies 
versus band splitting. As the equilibrium temperature of the 
system is O(eio), we can approximately neglect interactions 
and treat the final state as a free Bose gas. 



ture T"^ in a manner analogous to the virial theorem, 
just as Miiller et al. initially proposed fld\. To simplify 
matters, we consider that the trap is initially loaded with 
bosons just up to the point where double occupancies are 
favorable. Here the condition on the number of particles 
is given by = 2y^Uoo/A, where t/oo is the on-site en- 
ergy in the ground band and the trap potential is given 
by A{2x/X)'^. As the length scale of the trap is so much 
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longer than that of the optical lattice parameter we as- 
sume we can treat the dispersion semiclassically: 



(V.l) 



We treat interactions in the mean field. As the tempera- 
ture T'^ is much larger than the width of relevant bands, 
we ignore any momentum dependence of the interaction 
Unn' between particles in bands n and n'. Introducing 
the band relative occupancy Pn — jjJ2kfnk, the equi- 
librium mean field shift F^q in band n is 



^nO — PO X/ PmUnmPm, 



(V.2a) 



Pn 



4t/oo Efc {cxp [{enk + Fn-fi) - 1}'^ 
Ek Lii/2 {exp [- (e„fc + F„ - fi) /T-]}' 

(V.2b) 



with Life the polylogarithm. The Pn are determined self- 
consistently. Using these approximations, the equations 
for number, energy and relative occupancy p„ become 



= 1, 

n 

-E 

1 v ^ f°° dx 



dx 



N^J-oo A/2 



exp 



T 



N 



nk 



£nk{x) 



A/2 



exp 



= £1 



(V.3a) 
(V.3b) 

3 ■ 

(V.3c) 



The values for the equilibrium temperature T'^ and chem- 
ical potential fi obtained from these implicit equations 
are plotted in Fig. [181 The relative occupancies po-5 are 
compared to experimental results from (16] in Fig. 1191 
some differences between our theory and experimentally 
measured values are expected to come from neglecting 
excitations in the y— and z— directions. 



VI. LONG TIME DECAY 

We want to model the long time relaxation to the 
equilibrium computed in the previous section through 
a Boltzmann equation formalism. The full Boltzmann 
equation includes both spatial terms which lead to a 
spreading in the trap and momentum terms which equi- 
librate the different bands. To make our analysis math- 
ematically tractable we ignore the spatial terms, even 
though we suspect that they are important in under- 
standing experimental results. We treat the density 
within the trap as fixed by its equilibrium value and al- 
low thermalization only through momentum and inter- 
band relaxation. Rather than a complete theory of the 
later stages of evolution, this section should be consid- 
ered a first attempt to model the dynamics of interband 
scattering. 
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FIG. 18: (Color online.) The solution for the chemical poten- 
tial jj, and equilibrium temperature T'' including mean field 
shifts from Eq. (|V.3|) . Here we have plotted —fi/eio in order 
to easily compare the two graphs. Note the y-axis starts from 
0.5. 
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FIG. 19; (Color online.) The full probability distribution in- 
cluding mean field shifts from Eq. HV.3|I . The "X" markers 
represent the experimental measurement of the equilibrium 
values of bands n = 0, 1,2 [3]. The agreement with exper- 
iment is reasonably good as expected from an equilibrium 
calculation with proper Bose statistics; possible differences 
between our experiment and theory could arise from neglect- 
ing excitations in the y— and z— directions. 



Initially the trap is loaded just to the point where 
double occupancies become energetically favorable. As 
discussed in section |Vl this maximally loaded condition 
implies N = 2y/Uoo/A. We derive the Boltzmann Equa- 
tion first assuming classical statistics in Section FVI Al and 
then for Bose statistics in Appendix [F) We solve these 
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two equations numerically in Section IVIBI 



A. Derivation of Boltzmann Equation with 
Classical Statistics 

For classical statistics to be valid we require the occu- 
pation of each band to be relatively small. This condition 
fails for the lowest bands that have relatively high occu- 
pation. We relax this condition in Appendix |F] when we 
derive the Boltzmann equation with Bose statistics. 

We work in the semiclassical approximation where we 
can specify the occupation fnkj of band n, quasimo- 
mentum k, and site j. Although this seems at odds 
with the uncertainty principle, it is a valid approxima- 
tion when the wavelength of the optical lattice A/2 is 
much shorter than the band-dependent harmonic oscilla- 
tor length ^n, ~ (HX)^^^ (^Am';^) , with m';^ the char- 
acteristic mass scale in band n defined by the bandwidth: 
h'^/m'^X^ ^ tn- For the case of maximal loading, the ra- 
tio ^,i/A is roughly 



X 



A 



1/4 



(VI.l) 



For experimental parameters, C„/A ^ 1 for all bands 
but n = 0. For deeper lattice well depths, S,o/X> 1 sug- 
gesting the semiclassical approximation may not be valid 
experimentally for the lowest band. In the following, we 
make the simplification that the trap is sufhciently weak 
to apply the semiclassical approximation in every band. 
The semiclassical energy £nkj of a particle in band n, 
quasimomentum k and site j including the mean field 

n' k' j 



shift Fjij — ^ ^n'fc' ^nn' fv'k'-i IS 



^nkj 



£nk 



Aj" + F^^ 



(VI.2) 



Since the Fnj are small compared to the equilibrium tem- 
perature r*^, we can ignore the j— dependence and use 
Fnj ~ fno- The classical distribution of jnkj at thermal 
equilibrium is given by: 



nkj 



N 



exp{-enkj/F 



^E„fcjexp(-e„fej/T'=)^ 



(VI.3) 



where the equilibrium temperature T*^ is determined by 
conservation of energy. For the case of a maximally 
loaded trap this is given by: 



j\j 5Z ^"'^J -^"fci 

nkj 



N[e, + ^], (VIA) 



which simplifies to 

^ Uoo 

Tl , T.nk i^nk + F„^o) exp [- {Cnk + F^^) /T^] 

2 



exp [- (e„fc + F-^) /T-] 



(VI.5) 



where F^^q is the equilibrium mean field shift. The ^ in 
Eq. (|VI.5p is the virial term which represents the energy 
contribution from the trap potential. This term was first 
suggested to explain experimental data Neglecting 
all spatial derivatives, the Boltzmann equation becomes 



^fnikij 

dt 



ATS 

'^2 ns>n4 k2,k3,k4 



{fnikij fn2k23 fnsk^jfji^k^j) 5 

(VI.6) 



with the function ^nlntj given by: 

rXK^i' k2,h,ki) = ^ \u:^Z^ik,,k,, k,,k,)\^ X 

5 {ki +k2-k3~ki + G)x 

^ (Cnifeij + ^n2k2j ~ ^naksj ~ ^riikij) ■ 

(VI.7) 

We now make the simplifying ansatz that within a given 
band the occupation fnko is Boltzmann distributed, but 
the relative occupation of bands p„ can fluctuate: 



fnkO ~ PoPi 



exp [- (£„fc -I- fno) /T] 

A 



Pa = 



00 



(VI.8a) 



(VI.8b) 



where pq is the density near the center of the trap and T 
is an effective temperature determined by conservation. 
At equilibrium, the band occupancy probabilities Pn and 
mean field shifts Fno are given by: 



(VI.9a) 
(VI.9b) 



where the partition function Z = J^n ^nC^^)- We allow 
the temperature to fluctuate in time in order to enforce 
conservation of energy. Using the ansatz in Eq. (jVI.8[) 
and summing Eq. (jVI.6|) over ki gives an equation for 
the Pn's: 



dpm 
dt 



~Po E ^ntTAT) 

712 ,713 >n4 



PniPn2 



Pn^Prii 



A A A A 

(VI. 10) 



with ^nlnt (^) given by the expression: 



kik2ks /C4 



exp 



fnifci + -fniO ~^ ^n2k2 + ^n20 



(VI. 11) 
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Since the temperature T is much larger than the widths 
of the bands we are considering, the exponential factor in 
Eq. (|VI.11|1 depends very weakly on momentum. For a 
given well depth, the coefficients 7^^"'' (T) are nonzero for 
the values of n, 712, na, 714 where the two-body scattering 
process n+n2 ^3+714 conserves energy; see Eq. (|VI.7p . 
When the well depth becomes sufficiently large that the 



scattering process n + ^2 ^ 773 + 714 no longer conserves 
energy, 7"^"* will jump from a nonzero value to zero. 
These jumps have a substantial effect on lifetimes in the 
excited band, as can be seen in Fig. [201 

Expanding the Boltzmann equation [Eq. (|Vf .lOp ] to 
first order in deviations in both p„ and T from equilib- 
rium gives: 



dt 



Z ^ 

n2,n3>ni 



inin2 



5pm , 5pn2 Spn3 Spn 



A A 



-5T 



u„ 



u„ 



7140 



(VI.12) 



where () indicates the average over a single band 
J2nk 9nk exp {-enk/T) 



(gn) 



J2nk (-Enfe/T) 



(VI. 13) 



Since Eq. (jVI.lOp does not mix fnkj on different sites we 
must have energy conservation on-site. This requirement 
is expressed by the relation 



^ ^ (''"fc + ^no) fnkO — 0. 



(VI. 14) 



nk 



Eq. (|VI.14p leads to the following equation for 5T in 
terms of the 5pn '■ 



6T^ 



-E„5Pn((e„) + 2i^„-o) 



(VI. 15) 



We find the eigenmodes of Eq. (|VI.12p in Section lyTBl 



due to neglecting the spatial dependence of the Boltz- 
mann equation, which is known to be important in exper- 
iment but was neglected in our analysis. Inclusion of Bose 
factors gives only a small effect due to the high effective 
temperature T and relatively low band occupation prob- 
abilities Pn- The lifetime of the slowest-decaying nonzero 
mode of the Boltzmann equation with Bose statistics is 
shown for comparison in Fig. [20l We plot these lifetimes 
in units of the hopping time in Fig. [5TJ 

In Tableiniwe decompose the slowest-decaying mode vi 
of the linearized Boltzmann equation with Bose statistics 
onto fluctuations in the different bands: vi — ^„ vinSpn- 
The decomposition of the slowest-decaying mode of the 
linearized Boltzmann equation with classical statistics is 
nearly identical. For lattice well depths near IQEu the 
slowest-decaying mode has little overlap with Spi, sug- 
gesting that fluctuations from equilibrium in the first ex- 
cited band relax faster than fluctuations in other bands 
at this value of well depth. 



B. Solution of Linearized Boltzmann Equations 

We computed the {7^^^"" } with both classical and Bose 
statistics for nine bands over a range of well depths 
- AOEh; see Eqs. (|VI.11|) and (|F1T|) . Using these 
parameters, we determined the eigenvalues of the lin- 
earized Boltzmann equation for the cases of both clas- 
sical [Eq. (|VI.12P ] and Bose statistics [Eq. (|F.12j) ]: these 
eigenvalues represent the decay rates for each mode of 
the system. The linearized Boltzmann equation has two 
modes with zero eigenvalues corresponding to number 
and energy conservation. 

The lifetimes (1/rate) of the four slowest-decaying 
nonzero modes of the Boltzmann equation with classical 
statistics are plotted in Fig. [201 The comparison be- 
tween our numerical calculations and experimental mea- 
surements is relatively weak; this discrepancy could be 



VII. CONCLUSIONS 

We have considered the preparation, short-time evo- 
lution and relaxation of the p-orbital insulator state. A 
TT-pulse with characteristic time r sufficiently long that 
Ut ^ 1 will end in state which is a close approximation 
to the p-orbital insulator, with a density of excitations 
above the ground state which is exponentially small for 
large r. At short times, the insulator state is metastable 
due to the anharmonicity of the optical lattice. For well 
depths below the two-body threshold value Vt, our the- 
ory for the short time decay rate is broadly in agreement 
with measurements performed by Miiller et al. [4I ; above 
threshold there is a marked disagreement, perhaps due to 
a small fraction of double occupied sites from an imper- 
fect TT-pulse in experiments. Reducing this fraction could 
lead to substantially longer excited band lifetimes. As 



17 



-Classicall 

Classical2 
- Classical3 

ClassicaW 

Expt 

Bose1 



1+3^0+4 



1+2^0+3 




1+1^0+2 



Well 



20 30 
Depth (units of E ) 



40 











Vl2 


fis 


Vi4 


^^15 




5.0 




0.624 


-0.772 


-0.034 


0.110 


0.031 


0.005 


0.000 


10.0 




0.601 


-0.099 


-0.775 


-0.005 


0.156 


0.063 


0.013 


15.0 




0.452 


0.289 


-0.727 


-0.384 


0.118 


0.142 


0.056 


20.0 




-0.478 


0.854 


-0.145 


-0.129 


-0.065 


-0.013 


-0.001 


25.0 




-0.488 


0.850 


-0.166 


-0.098 


-0.043 


-0.021 


-0.008 


30.0 




-0.530 


0.832 


-0.014 


-0.142 


-0.078 


-0.027 


-0.011 


35.0 




-0.532 


0.833 


-0.029 


-0.130 


-0.071 


-0.027 


-0.013 


40.0 




-0.545 


0.826 


-0.011 


-0.109 


-0.081 


-0.029 


-0.016 



TABLE II: Eigenvector decomposition of the slowest-decaying 
mode of the linearized Boltzmann equation with Bose statis- 
tics [Eq. (|F.12|| ] on to bands n = — 6 for various well depths. 
The values for the case of a linearized Boltzmann equation 
with classical statistics are nearly identical. This table shows 
the evolution of the slowest eigenvector as various scattering 
rates change with lattice depth. For purposes of compari- 
son, the normalization condition Vi„ = 1 and a consistent 
sign convention were used. Note that vu passes through zero 
slightly above a lattice depth of WEh. 



FIG. 20: (Color online.) Comparison of the slowest-decaying 
modes of the linearized Boltzmann equation with both clas- 
sical statistics [Classical 1-4] and Bose statistics [Bosel] with 
the slowest-decaying mode extracted from published measure- 
ments of the occupancy of band n = 1 by Miiller et al. [Expt.] 
see Eqs. (|VI.12|l and (|rT2)) . The inclusion of Bose factors 
gives only a small effect, so we have plotted only the slowest- 
decaying mode for that case. As the optical lattice potential 
becomes deeper, some decay channels become energetically 
forbidden and their decay rates jump to zero. As noted in the 
figure, crossing some of these thresholds is accompanied by a 
substantial increase in the lifetime of the excited band. 
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FIG. 21: (Color online.) Lifetimes of the slowest-decaying 
modes of the linearized Boltzmann equation with both classi- 
cal statistics [Classicall-4] and Bose statistics [Bosel] in units 
of the hopping time; see Eqs. (|VI.12|I and (|F.12|I . This di- 
mensionless ratio is one of the most important figures of merit 
in characterizing the feasibility of studying quasi-equilibrium 
models in the first excited band. 



discussed in section IIIIl a single-tone Gaussian 7r-pulse 
may provide better fidelity than current experimental 
pulse shapes. 

The interband transitions lead to evolution of the ex- 
cited band occupation. This in turn affects the transition 
rates. The intermediate regime of relaxation is further 
complicated by a spreading in the magnetic trap, and the 
dynamics of this intermediate regime remains beyond our 
theory. At long times the system relaxes to an equilib- 
rium gas of nearly free bosons, where the chemical poten- 
tial and temperature describing the gas are determined 
by energy and number conservation. We modeled this 
asymptotic relaxation using a Boltzmann formalism with 
both classical and Bose statistics. The calculated decay 
rates differ considerably from experimental values; this 
disagreement could come from neglecting spatial varia- 
tions in the trap. 

The double-well optical lattice could offer an improve- 
ment for excited band lifetimes in the insulating regime. 
The double well lattice can satisfy the hard anharmonic- 
ity condition, £2 + ~ 2ei > 0, for suitable lattice pa- 
rameters. This condition changes the initial decay mech- 
anism from an exothermic process to an endothermic one, 
leading to a significant enhancement of the excited band 
lifetime. Stojanovic et al.'s work in the superfluid regime 
found that excited band lifetimes in the double-well lat- 
tice could be orders of magnitude larger than equivalent 
single- well lattices [l3| . 

Our theory provides a framework to understand a wide 
range of experiments performed in the excited band, in- 
cluding those performed by Miiller et al. We hope 
these insights and calculations can improve the experi- 
mental realization of excited band models in optical lat- 
tices. 
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APPENDIX A: SCHRODINGER EQUATION OF 
TT-PULSE IN TWO-SITE MODEL 



As discussed in Section lllli we can understand many 
aspects of the preparation of the p-orbital insulator using 
a model of two bosons in a well of two sites. Each site 
has an s level and a p level, but we only include tunnel- 
ing in the p-band. We initially prepare the system in the 
two-particle ground state, and then promote both parti- 
cles to the excited level through a 7r-pulse. The model 
Hamiltonian H is given by : 



+ -y^^LO^LO^iO^'LO + Uwb\^b\^bLibLQ 
+ {L^R). 



(A.l) 



The six relevant parity-symmetric eigenstates are given 
to 0{ti/U) in Eq. pil.2[) . We promote the bosons from 
the s-band to the p-band using a pulse of the form: 



pulse = A(i) e-'"^' (bl^bLo + b^RibRo 



h.c, (A.2) 



with the Gaussian envelope A(t) = ^ exp(— t^/r^). We 
choose the dimensionless pulse strength a to optimize a 
single TT-pulse . The excitation can be considered a stan- 
dard two-photon process, with the frequency lo chosen to 
be half the energy difference of states |1) and |5). 

The Hilbert space naturally divides between a sub- 
space composed of states with nearly one particle per site, 
{|1) , |3) , |5)} and a subspace composed of states which 
include doubly occupied sites {|2) , |4) , |6)}. To under- 
stand the behavior of the Schrodinger equation, we first 
solve for the evolution within the {|1) , |3) , |5)} subspace 
and then use that time evolution as a drive term for the 
{|2) , |4) , |6)} subspace. 



1. Evolution of {|1) , |3) , |5)} subspace 

For the following we denote the amplitudes of the 
|1) , |3) , |5)} by a = {ai, 03, 05}. The Schrodinger equa- 
tion in the [/ — > 00 limit becomes 



au 



(A.3) 



with Hamiltonian -ffi35 given in the frame rotating at 
frequency lo by: 











\/2ae-"' 




\/2ae 





\/2ae-"' 






(A.4) 



We operate in the regime k <C 1, where we can expand 
the solution a as a power series in k; see Eq. (jIII.6|) for a 
definition. Changing variables from u to v: 



(A.5) 



V ^ a^j -{1 + cri u) , 



where erf denotes the error function. Here v ranges be- 
tween and a\/2T:. In terms of these new variables, 
Eq. (|A.3|) becomes 



dv 







1 











;)•( 




^0 


1 































a, (A.6) 



where u is an implicit function of v. Since Eq. (jA.6p 
obeys the constraint 01—05 = 1, we can simplify the 
Schrodinger equation in terms of the new variable z: 



dv 



2 
2 

-\/2 
\/2 



1 -1 
-1 1 



(A.7a) 
(A.7b) 



Expanding z as a power series in k, we find both the first- 
order corrections /i and /2 as well as the second-order 
corrections hi and h2'. 







h{v) 


~ 2 V2a io ' 


f2{v) 




hi{v) 






e"'(''")cos 






h2{v) 


=K{v). 



(A.Sa) 



(A.8c) 



V2{v' - v") sin (^V2 v" ] , (A.Sd) 

(A.8e) 
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In terms of the functions fi{v) and hi{v), we can express 
{01,03,05} to 0(k2): 



ai{v) wcos^ (\^) + «Im[/i(u)]K + Re[/ii(w)]K^, 



(A.9a) 



asiv) « - V2Re[h{v)]K 

a^iv) « - sin^ ("Tf ) ^ + Re[/ii(w)]K^ 



(A.9b) 
(A.9c) 



2. Calculating 05] after application of pulse 

We now determine the fraction of particles excited to 
state |5), the analogue of the p-orbital insulator, in the 
limit K <C 1. We use the transformed variable v defined in 
Eq. (jA.Sp . For k = 0, |a5(Q;\/27r)| reaches its maximum 

value of unity when a = Expanding a as a Taylor 
series in n gives 



a = — — \-bK, + CK + 



(A.IO) 



In terms of these coefficients, the leading order correction 
to lasl^ after the pulse is 



a5(Q;v 27r) 



05 



V2 



2Trb^K'^ + 0{k''). (A.U) 



To 0(k ), |a5(Q;v27r)| is maximized when b — 0. Our 
final expression for the occupation jasl^ after the pulse is 



05 



V2 



= 1 - 



fli 



V2 



03 



V2 



1 - (0.265... )k^ + 0(k^ 



, (A.12) 
(A.13) 



3. Calculating the dimensionless pulse strength a 

In order to determine the dimensionless pulse strength 

a to 0{k^), we need to maximize |a5(a'\/27r)| over the 
coefficient c in Eq. (jA.lOP : 



a5iaV2^) =1 - (0.265... 



(^cV2Tr + f3^ - 7c 



where the values (3 and 7 are given by: 

/3 = 0.0825... , 
7 = 0.9476... . 



(A. 14) 

(A.15a) 
(A.15b) 



To 0(k ), a5(av27r) is maximized when: 



1^^ = 0.0425, 

4tt 



(A.16) 



The optimal value of the dimensionless strength a which 
maximizes the occupation of state |5) is given to 0{k^) 
by: 



a = - (0.0425 . . .)k^ + O {k^) , 
The corresponding value of jasp is 



(A.17) 



a5{aV2^) = 1 -(0. 265... + (0.0045... )k'* + 0(k^). 

(A.18) 

4. Evolution of {|2) , |4) , |6>} subspace 

To estimate the amplitudes 02, 04, and ag of the 
{|2) , |4) , |6)} subspace, we substitute the leading terms 
for amplitudes {01,03,05} from Eq. (|A.9|) into the 
Schrodinger equation. This procedure is valid to lead- 
ing order in the small parameter k/U: 



da2 ... , ^/7^^l 

— — « -1(7202 H 1= ' 

du 2V2?7io 



C?04 

du 



'V2U- 



' e-"' cos^ 



cos (■^c'l'f J (A. 19a) 



— (1 + erf u) 
4 



2TTt 



1 / 1 ^11 \ -u-' ■ 2 
1 — Trrr- e sm 



(1 + erf u) 



du 



Un V 2t/io, 

(A.19b) 

(A.19c) 



with {U2,U4,Ue} defined by 



Un 
Ue =C7 + 



8k 



r, Uii ■ 



(A.20a) 
(A.20b) 
(A.20c) 



From Eq. (|III.6p . we make the connection between ti/Uu 
and k/U : 



ti 

Un 



We expand {02, 04, og} to leading order in k/U: 



(A.21) 



|a„(oo)|" 



O 



k/U 



if ra = 2,4,6 . 

(A.22) 
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Solving Eq. (jA.19p gives us the expressions for with r = Uh/Uiq. In the limit of large C/, 03"^ \ 04"^' and 
{aW,a«,a«}: 



;(2-r) 



due 



-u +iU2U 



COS 



/ Trerf 



(A.23a) 





2 27r 


/ du< - cos 




^ 2-r 


i-oo 1 2 



TT (1 + erf u) 



,(1) 



1 - -J sin^ 



TT (1 + erf 



-u -\-iU4u 



(A.23b) 



/ Trerf ^ 



V 2 



(A.23c) 



flg^'' are all exponentially suppressed. We can approxi- 
mate these terms using the stationary phase approxima- 
tion (see Appendix [B|) : 



102(00)!' 



|a4(oo)|' 



|a6(oo)|' 



8U i2Ulo - UnUw) 2:0(1/2) 
ttkUio C/4 exp [— 2[/4a;o(J74)] 



U2 exp [-2[/2Xo(f/2)] ^ ( U2 

{2xo{U2) 



2t/(2C/io-i7ii) 



Xo{U4) 



2x0 (C/2) 



C/io/ l2a:o(C/4) 



U{2Uio-Un) 



mo J 



-erfi xo{Ue 



2xo{Ui 



r 



(A.24a) 
(A.24b) 
(A.24c) 



with the fimction xo{U) given by: 



xo{Un) = Wlog-^. 



(A.25) 



In Fig. [55] we compare numerical determinations of 
hnif^^o \an\^ / K with a stationary phase approximation 
calculated in Eq. (jA.24p Since Uq < U2,U4, og 02,04 
due to the exponential dependence on [/„. This makes 
Og the most important error in determining |a5(oo)| : 



105(00)!' 



K 
U 



,(1) 



(0.265 . . .)k^ + (0.0045 . . + 0[k 



(A.26) 



APPENDIX B: STATIONARY PHASE 
INTEGRATION 

The following integral frequently arises in our discus- 
sion of the preparation of the p-orbital insulator using 
TT-pulses with Gaussian envelopes: 



Y{k) 



dx exp 



ikx — i— erf(a;) 



(B.l) 



We want to derive an asymptotic expansion for Y{k) in 

the limit k ^ 1. We define the inverse function for y{x) = 

2 

e"^ at all points in the complex plane by choosing the 
branch cut to be along the negative imaginary axis: 



y ^{x) = i\/\ogx. 
Y{k) has a stationary phase point at xi: 



Xi 



(B.2) 



(B.3) 



This formula is valid for k both positive and negative. 
Using the standard stationary phase approximation, we 
can write an expansion for asymptotically large k: 



Y{k) 



Vk 



1 



exp 



!fc!Jiog— X 



21og^ 41og^^ 



(B.4) 



From Eq. (|B.4[) . we see Y(k) decreases faster for large 
negative k than large positive k because the magnitude 
of the logarithm term is larger: 



log- 



= Jlog^ 



!fc! 



(B.5) 
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FIG. 22: (Color online.) Comparison of the numerical calcu- 
lation of lim„^o \an(oo)'^ / K for n = 2, 4, 6 with the stationary 
phase approximation; see Eqs. (|A.23|I and (|A.24|I . 



A more qualitative way to understand this result is to 
expand the exponential of F(fc)'s integrand in Eq. (jB.ip 
for small x: 

— + ikx — i— erf(x) « —x'^ + i{k — \pn^x. (B.6) 



The integrand is more slowly oscillating for fc > than 
for fc < 0. 



1. Related Integrals 

In the course of computing the density of excitations 
in section UlIBl the following two integrals c(fc) and s(fc) 
arise: 



c(fc) 



due 



iku—u 



cos I —erf u 



Y{k) + Y{-ky 



due"'''-'' sin 



-Y{k) + Y{-ky 
2i 



gerfu) 



(B.7a) 

(B.7b) 
(B.7c) 

(B.Td) 



Both integrals share the same stationary phase point xi: 



(B., 



In the stationary phase approximation, the integrals be- 
come 



c(fc) 



5(fc) =i 




(B.9a) 



(B.9b) 



APPENDIX C: DERIVING RECURSION 
RELATIONS FOR Pg" AND pexc 



In Section llll B( we consider the problem of loading the 
excited band in an system of N sites. Naive perturbation 
techniques fail because the expansion parameter changes 
from ti/U to ti^/N/U, a large parameter in the ther- 
modynamic limit. We get around this difhculty by first 
allowing hopping only between n sites, where n is suffi- 
ciently small that ti^/n/U <C 1 and perturbation theory 
remains valid. We derive recursion relations in the vari- 
able n for two important quantities: the probability Pg 
that the final state after application of the pulse is the 
p-orbital insulator and the "density of excitations" poxc 
of the final state defined in Eq. (|III.19p . By solving these 
relations for n N, we derive expressions for the prob- 
ability Pg and the density of excitations pcxc accurate in 
the thermodynamic limit. 

An excitation pulse A(t)e*'^* promotes bosons from the 
ground band to the first excited band (see Fig. 2]). We 
choose the pulse carrier frequency uj to maximize the oc- 
cupation of the p-orbital insulator. We are interested in 
pulses sufficiently long that U = Ut ^ 1, so that the ex- 
citation density above the ground state is small. In the 
following discussion, the phase angle e{t) = dt'A{t') 
frequently arises. 

By changing to the interaction picture, we transform 
the time dependence of the pulse into the time depen- 
dence of the tunneling ti and the pulse detuning eig — lo. 
Treating all of the Hubbard U parameters as equal, the 
Hamiltonian with tunneling only allowed between the 
first n nearest neighbor pairs [see Eq. (jIII.20[) ] is given by: 

n 

n 

-t- (eio - uj)'^{c b\,j + is [c bij - is b^j) + h.c. 



U 



(C.l) 
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where c = cos 9 and s = sin 6*. Initially, we prepare the 
system in the true iV-particle ground state |?/'o) at a den- 
sity of one particle per site: 



N 



(C.2) 



The wave function j^*") describes the evolution of the 
system under i7„. Expanding to leading order in 
the small parameters ti^/n/U and tiTy/n/U gives 



N 



1*7(0) =A4t) Yl bl^, |0) 



n TV 



N 



nt] ntfr 



1 '"-i' 



C/2' C/2 



(C.3) 



The coefficients B^^\ C'"'^' and must all be indepen- 
dent of n by translational invariance, whereas An must 
decrease with increasing n in order to preserve normal- 
ization. Expanding the Schrodinger equation for B^^\ 
C*-^-*, and I?^^^ to leading order in small parameters gives 

d(i) 



— « [7 fiV2 sin 6* cos 61 + 5(1)) , (C.4a) 
du \ J 

(C.4b) 



^AcW «;7f2sin2 + C(i 
du V 



sin6'cos6' + \/2cos2 ^s'^^ 



z sm y cos - z ( ^^^^ ) cos^ BD^^'^ . 
tj/U 



(C.4c) 



Multiplying both sides of Eqs. ((a4a| and ([Olb)) by e*^" 
and integrating by parts gives 



b'-^'^ = - iV2sm9cos9 + 5B^^\ 
=-2sin2 + <5C7«, 



(C.5a) 
(C.5b) 



where the terms (5_B'^^^ and SC'^'^^ are given by: 

du'^ cos 26*6*^"', (C.6a) 

-OO 

5&^^ =2e-^^'" ! du'^sin26le^^"'. (C.6b) 



The terms 5B'^^^ and 5C^^^ in Eq. (|C.5|) represent exci- 
tations above the ground state and so contribute to the 
density of excitations Poxc- Substituting our solutions for 
5(1) and C(i) from Eq. ([aS]) into Eq. (jOl)) . we get the 
following expression for £)(i)(oo): 



^ - eio 



h 
h 



du s\nO{u) cos^(w), 



(C.7a) 
(C.7b) 



[2 cos^ 9{u) cos 29{u') + sin 26'(u) sin 29{u')\ . 

(C.7c) 



To minimize |£''-i)(oo)| , we choose the frequency uj to 
cancel the real part of D^^'>{oo): 



Hi 

U 



Rejh 
h 



(C.8) 



For [/ 1, Eq. (|C.8|) reproduces the two-site result of 
Lo Ki tiQ — 2t\T jU . For this optimal frequency, the value 
of £'(i-'(oo) becomes 



U{u'-u) 9{u')x 





2 


/"OO I^U 




/ du du' sin 









{cos26I(m') +COS [26I(m) - 26I(m')]} 



(C.9a) 



We have plotted [D'^^) (cx)) | for the specific case of a 
Gaussian pulse in Fig. [23l 

We can determine the ground state wave function | 
of Hn{oo) by enforcing the adiabatic pulse condition: 



1*;^) = Um m) 

U^oo 
N 



(C.lOa) 



N n N 



f2 \ 



(C.lOb) 



In order to compute the probability Pg of completing the 
pulse in the p-orbital insulator state and the density of 
excitations Pcxc, it is convenient to reexpress |^'"(oo)) 
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from Eq. (|C31) using |*^): 



0.4r 



u 



n 



(C.ll) 



To derive the recursion relations, we first calculate the 
change in An as n — > n + 1 to preserve normalization: 



An+l{oo) 



1 



2C/2 



^2 

2C/2 



C/2 



4 +6^2 



(C.12) 



Taking the modulus square of Eq. (|C.12|) gives us the 



recursion relation for P" 



pn+l 

log — 



''1 

24t2 

c/2 



o 



C/4' C/4 



(C.13) 



To derive the recursion relation for the density of excita- 
tions Pexc) we use the definition in Eq. (jlll.lQp and the 
wave function 1^-7+^) from Eq. ([aTT|) : 



Pcxc{n + 1) =Pcxc(n) + SB^^'> {00) 



O 



/ f4 
( 1 



Vivc/4 



(C.14) 



Note that we do not include |£'(^)(oo)|2 in the expression 
above as it is 0{k^), negligible in experimentally inter- 
esting regimes. 



APPENDIX D: MODEL DEMONSTRATING 
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FIG. 23: The coefficient ]d(i)(oo)] for the specific case of 

a Gaussian pulse; see Eq. (IC.7|I . As expected, it is strongly 
suppressed for large Ut. 



state \i'Bf)- We denote the decay rate to \i'As) by Ri 
and the decay rate to j^'B/) by i?2- In general, R2 is 
0(ii/C/)2 smaller than and so can be neglected. We 
demonstrate this with a simple three-site model that has 
the essential features of the full N-site system. 

Our model consists of three sites at j = —1, 0, 1, each 
with n = and n = 1 band levels. The sites are linked 
by nearest neighbor tunneling ti in the excited band. 
The on-site interaction between particles in bands ni and 
712 is Un-in^ In addition, we treat the n = 2 band as a 
wide continuum, indexed by the quasimomentum q. To 
understand the 1-1-1 0-1-2 decay channel, we introduce 
an interband term Uil that scatters a pair of central site 
n = 1 band particles into an n = 2 band continuum state 
and a localized n = band particle. The Hamiltonian 
neatly divides into two pieces: -ffo which preserves band 
index and Hi which leads to interband transitions. In 
second quantized notation, and Hi are given by: 



J=0,±1 q 



bii 



h.c 



+ h [bio {b 
fb\,h\ 

Hi=U"i'i-^Y.^lbl {biof+h.c, 



-^^10^10^10^10 + C/io^oo&io ^00^*10, (D.la) 



(D.lb) 



In section IIV| we discussed that the p-orbital insula- 
tor can directly decay via the 1 -f 1 0-1-2 channel 
into two energetically permissible final states, \tpAf) and 
I'f/'B/) distinguished only by the nature of the holes in 
the n = 1 band. One n = 1 hole is localized around the 
n — particle in a bound state in state IV'A/); whereas 
the n = 1 hole and n = particle are well separated in 



where the operator 6^^^ creates a particle in band n at 
site j. We initially prepare the system in \ipo), the lowest 
energy eigenstate of /Cq that contains three particles all 
with band n = I indices. This eigenstate is the model 
analog of the p-orbital insulator. To leading order in 
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ti/U, the eigenstate lipo) and its energy are given by: 



l^o) 



4t2 



J7- 



11 



|0), 
(D.2a) 

(D.2b) 



We denote the Hq eigenstate analogous to ji/'A/) by \'i{}2q) 
and the eigenstate analogous to \tpBf) by IV's?')- To lead- 
ing order in ti/U, these eigenstatcs and their energies are 
given by: 



fell + 61-1 
V2 



hV 2 

Uu 



biib 



1-1 



61161-1 + 77-610 (611 + 6i_i) 
1^11 



•2q} 



£20 



ei + Co 



2t|_ 

C/10^ 



(V'Sg'l -f^O IV'Sg') ~ £29' + ei + eo + Uio + 



2t2 



10 



(D.3a) 

iV-o), 
(D.3b) 

(D.3c) 
(D.3d) 



The width 4t2 of the n = 2 band is sufficiently large 
that both \'4'2q) and IV'Sq') conserve energy for specific 
values of 9,9'. We use Fermi's Golden Rule to compute 
the decay rates Ri and i?2 from to \'4>2q) and IV's?') 
respectively. 

Of 

(V'2,|i?i|V'o)«Tr^t/°f, 



(V^s^'lifilV-o) 

|(^2g|i?l|V'0)| 



t/ll 

-2V2t^ 



2 



rr02 

^^11 7 



-Ri 



i?2 



(D.4a) 
(D.4b) 
{D.4c) 
(D.4d) 



Since the matrix element {tpsqi \Hi\^o) is 0{ti/U) 
smaller than {ip2q\Hi\ipo) , the rate R2 is 0(ii/C/)^ 
smaller than the rate Ri and can be neglected. This 
insight holds for the full A^-site problem considered in 
Section HVl 



1. Two-body decay rate as V ^ 

To evaluate the two-body decay rate near threshold, 
we use the tight-binding approximation with mean field 
shifts and neglect hopping in band n = 0. The band 
dispersion in the lowest three bands is given by: 

eokiV) « eoiV), (E.fa) 
eikiV) « eiiV) + 2hiV) cos(7rfci), (E.fb) 
hkAV) ~ hiV) ~ 2<2(V^)cos(7rfc2). (E.fc) 

The Fermi's Golden Rule formula for the two-body short 
time decay rate F is 



F 



X S (e2 



dki cos Trfci / dk2 Uif {k 



2t2 cos 7rfc2 



2ei — 2ti cosTrfci) 



(E.2) 



We denote the maximum of the delta function argument 
&sg{V): 



9{V) = €2 + 2i2 + eo-2ei + 2ti. 



(E.3) 



If g{V) > 0, there exists a region in parameter space 
where the delta function is satisfied and the two-body 
rate F is nonzero. For g{V) < 0, the delta function's 
argument is never zero and F vanishes. The condition 
<?(Vt) — sets the threshold value V*. For the experi- 
mental parameters in [3], Vt = 25.97 Eji. 

For well depths near threshold, the integrand is 
nonzero only for k2 ~ 1; as Uii{k) is relatively insen- 
sitive to fc, we can pull |/7°^(A:2 = 1)| outside the inte- 
gral. Evaluating the delta function and substituting the 
variable y — \/^tilg cos gives: 



where the prefactor Fq is given by: 



3/2 



ri/2- 



(E.4) 



(E.5) 



APPENDIX E: DECAY RATE AT SHORT TIMES 
NEAR TWO-BODY THRESHOLD Vt 



Since g{V) ^ ti close to threshold, we can expand the 
two-body decay rate as a power series in g/ti: 



We discuss the behavior of the decay rate at short times 
for well depths near the two-body threshold value Vt- 
We approximate this decay rate using the two-body rate 
1 -f 1 ^ -I- 2 for well depths below the threshold value 
Vt and the three-body rate f + f + f^0 + + 4for well 
depths above threshold. As the well depth V — > Vf' , 
the two-body rate rises linearly to a constant. However 
for V Vf^ , the three-body rate instead diverges as 
[iV-Vt)/h]-\ 



-Fr 



dy 



: X 



g_ 

ti 



8t2 



o 



2 2 

<r_ <r_ 

i? ti. 



(E.6) 



We treat the n — 1 band hopping energy ti as much 
smaller than the n = 2 band hopping energy t2, consis- 
tent with experimental values. Taylor expanding g{V) 
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to leading order in {V — Vt) gives the two-body rate be- 
havior near threshold with the leading order correction 
in iV~Vt)/ti: 



T{V < Vt) =ro 

r{v > Vt) =0. 



1 



7 dg{V') 



16 dV 



V^Vt 



, (E.7) 
(E.8) 



2. Three-body decay rate as 1/ — » V^^ 



We now consider the three-body decay rate for well 
depths just above Vt- As V ^ Vt^ , the three-body rate 
diverges because a virtual intermediate n = 2 state be- 
comes progressively closer to a real, energetically permis- 
sible state. The three-body rate is given by: 



1287rt? ^ dki 



dki 



cos^ (Trfci) X 



1 dk U^l{k)U^I{k,ki) 



/-I 2 eifc -I- ei — £0 — e2fc 



(E.9) 
(E.IO) 



Near the singularity, the dominant contributions to the 
rate T come from momenta near (fc, fci) — (1, 1). Pulling 
nonsingular terms like Ui1{l) outside the integral and 
exploiting left-right symmetry gives 



j^Jo J -I 2 eife, +ei-e, 







1) 


2 


u'Ai^X) 


2 




di 


fe4/*4|fe0 





Co ^ e2/c 

(E.lla) 

(E.llb) 



where k1 is set by conservation of energy: 
£4^.0 -I- 2eo - ei - en = 0. 



(E.12) 



Using the parameter g{V) from Eq. (|E.3[) . the fc— integral 
in Eq. (|E.11[) can be performed using a contour: 



1 dk 



_i 2 eifei + ei - Co - £2*; 2t2WP' 
^ g-4tiCOS^(4i) 

2*2 



(E.13a) 
(E.13b) 



Using the standard approximation 1 — a « exp(— a) for 
a <C 1 and expanding — 1 to leading order in the small 
parameters g/t2 and ii/t2, Eq. (jE.lip becomes 



ri exp 



2W(-5 + 4iiCos2 2|i)/t2 



4*1 cos2 



(E.14) 



Summing the geometric series gives 



32^^^^ 



Trfc 



cos 



-3/2 



(E.15) 



If 5(y)/*i < 0, as it is for V >Vt, this integral exists and 
is finite. For V < Vt this integral diverges reflecting the 
presence of a real intermediate state. We can express the 
solution to Eq. (jE.15p using the complete elliptic integral 
of the second kind with imaginary argument: 



dx 



K + cos' 



_3/2 2 E (i/VK 



K 



, (E.16a) 



2 



O(logX). 

(E.16b) 



Substituting —g/Ati for if in Eq. (jE.16P we have the final 
form for the three-body rate F near the threshold value 
of well depth Vt: 



8(ti)^/' 


\u?Hi)\'\u?i{i,k2) 


2 




d€4k 

dk 


1,0 
K4 


d 
dV 




~Vt) 

( 



(E.17) 



APPENDIX F: DERIVATION OF BOLTZMANN 
EQUATION WITH BOSE STATISTICS 

In this Appendix, we derive the linearized Boltzmann 
equation at long times assuming Bose statistics. The 
principles of the derivation are nearly identical to those 
of Section fVl Al but the formulas are notably more com- 
plicated. We work in the limit of a weak trap, where 
the semiclassical approximation is valid. We specify the 
occupation fnkj of band n, quasimomentum k and site j 
with the corresponding energy enkj approximately given 
by: 



£nk + Aj' 
1 ^ 



nj I 



^nj — 7 , Unn'fn'k'j, 
n' k' 



(F.l) 
(F.2) 



where Fnj is the mean field shift. Since the Fnj are small 
compared to the equilibrium temperature T'^ we can ig- 
nore the j— dependence and use Fnj ~ F„o- The Bose 
distribution of fnkj at thermal equilibrium is given by: 



f:,^^ {exp [{enkj~f^)/n 



(F.3) 



where the chemical potential /i and equilibrium temper- 
ature T'^ are determined by simultaneously satisfying en- 
ergy and number conservation: 



^ f-nkjfnkj - ^ ( ^1 + ^ 



1 



nkj 



^J2fn,,=N. 

nkj 



(F.4a) 
(F.4b) 
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To simplify the discussion we introduce z^f. 
exp [- (e„fc + F„o - /T"]- Eq. (|R4a)) simplifies to 



Enfe (enfe + ^»o) Lil/2 Kk) + ¥^13/2 «i 



Z]«fc Lli/2 (z^j.) 



El 



fA)0 

3 ' 
(F.5) 

where Life is the polylogarithm; see Eq. (jVI.Sp for com- 
parison. Neglecting spatial derivatives, the Boltzmann 
equation with these approximations becomes: 



^fni kij 

di 



= -E E ^3 E KlZAh,h,k„k,)x 

"2 n3>n4 k2,k3,ki 

[fnikijfn2k2j (1 ~^ fnskaj) (1 + fn4,k4j} 
~ fnsksj fn4,k4j (1 + fnikij) (1 + .fn2k2j)] j 

(F.6) 



with the function T^-J^J^. given in Eq. (|VI.7p . We are most 
interested in scattering near the center of the trap, and 
so we set j — for the remainder of the discussion. We 
make the ansatz that within a given band the occupation 
fnko is Bose distributed, but the relative occupation of 
bands p„ can fluctuate. This ansatz is given by: 



where the parameter T is an effective temperature set by 
conservation of energy. At equilibrium, the band proba- 
bilities Pn and mean field shifts F„o are given by: 



FnO = E Unn'dniT"). 



(F.8a) 
(F.8b) 



/nfcO 



^ {exp [(e„fe + F„o - m) /T] - 1}"' , (F.7a) 

^T"" Ek Lii/2 {exp h (enfc + F„o - /i) /T^]} ' 

(F.7b) 



(F.7c) 



Using the ansatz in Eq. IF. 71 and summing Eq. (jF.6[) over 
ki gives an equation for the p„'s and T: 



dt 



1 - ^ 



E 

n2 ,n3>n4 k\k2 k^k4^ 



'^'mP'n2PniPn2 + d^ij ^'"a ^^"3 fcs) (l + P™4rfn4fc4) 



l + d 



■n^k^ 



PniPniPnsPrii 



l + d 



'n2k2 



(F.9) 



where d„fc is given by: 



its band average d„. We define ■^^f^t by 



dnk = {exp[{enk + F„o- fi)/T]-l} \ (F.IO) 



min2 



74 E ^«?"20 (fcl-4)rfnidn2 (l+d„3)(l + '^n4)- 

(F.ll) 



iV4 



fci- 



Since the temperature is so much larger than the band- 
width of the relevant bands, we can approximate dnk by 



Expanding the Boltzmann equation [Eq. (jVLlOp ] to first 
order in deviations in both p„ and T from equilibrium 
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including these approximations gives 



dt 



~l~ ^ ^ ^Pni (^mni ~l~ ^nin2 ^mn^ ^771714) 
m 

(F.12) 

where the terms {Vny^niBmn} are defined by 

_ 1 1 ^ + F'=„ - //"i 



^rnn 



T-df,il + df,)N 



(1 - 



(F.13b) 

ly^^ (F.13c) 



with and are shorthands for d„(T'^) and A„(T^) 
respectively. Since Eq. (|F.9|) does not mix fnkj on differ- 
ent sites we must have energy conservation on-site, giving 
the relation 

6Y,(^nk + F^o)fnko^0. (F.14) 

nk 



Eq. (|F.14p leads to the following equation for ST in terms 
of the Spn'. 



6T 



with (e„) given by: 



(F.15) 



(F.16) 



We find the slowest-decaying eigenmode of Eq. (|F.12|) in 
Section fVl Bl The decomposition of the slowest decaying 
mode Vi on each of the Spn is given in Table |TT1 
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